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A new mechanism for toroidal momentum transport in a tokamak is investigated using the gyro- 
kinetic model. First, an analytic model is developed through the use of the ballooning transform. 
The terms that generate the momentum transport are then connected with the poloidal derivative 
of the ballooning envelope, which are one order smaller in the normalised Larmor radius, compared 
with the derivative of the eikonal. The mechanism, therefore, does not introduce an inhomogeneity 
in the radial direction, in contrast with the effect of profile shearing. Numerical simulations of the 
linear ion temperature gradient mode with adiabatic electrons, retaining the finite p* effects in the 
ExB velocity, the drift, and the gyro-average, are presented. The momentum fiux is found to be 
linear in the normalised Larmor radius (p.) but is, nevertheless, generating a sizeable counter-current 
rotation. The total momentum flux scales linear with the aspect ratio of the considered magnetic 
surface, and increases with increasing magnetic shear, safety factor, and density and temperature 
gradients. 

PACS numbers; 52.25.Fi, 52.25.Xz, 52.30.Gz, 52.35. Qz, 52.55.Fa 



I. INTRODUCTION 

For a tokamak, plasma rotation can be beneficial for 
confinement and stability. In particular, if the rotation 
is associated with a suflSciently large ExB shear, it sta- 
bilises turbulence [H-Il], and leads to an enhanced en- 
ergy and particle confinement. For this reason, there is a 
strong interest in the experiments on, and the theoretical 
description of, plasma rotation. An observation of par- 
ticular relevance to a fusion reactor is the occurrence of 
spontaneous rotation i.e. a plasma rotation that 

develops vifithout an externally applied torque. The the- 
oretical description of this spontaneous rotation has at- 
tracted much interest, and significant progress has been 
made in recent years. For an overview we refer the reader 
to Ref. 

All contributions to toroidal momentum transport are 
associated with a parity symmetry breaking mechanism 
in the gyro-kinetic equation |15l4l7l |. It was shown in 
that to lowest order in the normalised Larmor ra- 
dius = pi/i?o ^ Ij where pi is the ion Larmor radius 
and Ro is the major radius of the magnetic axis, a finite 
momentum flux is generated only in the presence of a 
rotation gradient 1 181 or for a rotating plasma through 
the Coriolis pinch [19|, [l^l or for a non up-down symmet- 
ric magnetic equilibrium pTI . [22| . In next order many 
different mechanisms appear. The most studied mecha- 
nisms to day are: the effect of neo-classical equilibrium 
flows [1^, the ExB she ar |2'4n 28il and some form of ra- 
dial profile shearing |29l - [3l| . The latter two mechanisms 
introduce an inhomogeneity in the radial direction. 

In this paper a new mechanism of toroidal momentum 
transport is studied. The mechanism is connected with 
higher order, in the normalised Larmor radius p*, paral- 
lel derivatives of the perturbed distribution as well as the 



perturbed fields. For micro-instabilities, the perturbed 
quantities have scale lengths perpendicular to the mag- 
netic field that are of the order of the Larmor radius. The 
fast transport along the magnetic field, however, strongly 
damps short scale perturbations, and the scale length 
along the magnetic field is of the order of the system size 
(Ro). This ordering allows some parallel derivatives that 
appear in the model equations to be neglected since they 
are one order smaller in the normalised Larmor radius 
(p*) compared with the perpendicular gradients. The 
neglected terms, however, do break the parity symmetry 
and might, despite their smallness, be of relevance in the 
description of momentum transport. In this paper the 
effect of these parallel derivatives on momentum trans- 
port is studied. The paper concentrates on quasi-linear 
theory, leaving the non-linear state for further study. 

Before turning to the complete model description and 
its numerical solution, the physics mechanism is dis- 
cussed below using an analytic model. We stress that 
this model does not include all contributions to the mo- 
mentum flux obtained in the model that is developed in 
the next section. The goal of the analytic model is to 
clarify the physics and to provide a first estimation of 
the expected momentum transport. The physics mecha- 
nism considered in this paper can be most easily clarified 
using the ballooning transform with s — a geometry [3^ . 
All perturbed quantities are then assumed to have the 
form 

G(r, 0, if) = G{e) exp[in(v3 - qir)0) - iwt] + c.c, (1) 

where ip [9) is the toroidal (poloidal) angle, r is the radius 
of the flux surface, n 3> 1 is the toroidal mode number of 
the instability, q is the safety factor, lj = lob. + 17 is the 
complex frequency, and c.c. denotes the complex conju- 
gate. The rapid variation (n ^ 1) of the perturbation 
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perpendicular to the magnetic field is represented by the 
eikonal, with the argument ((^ — qO) chosen such that it 
is constant along the magnetic field. The envelope G is 
then assumed a slowly varying function of the poloidal 
angle {9). The ExB velocity (v^ = b x V(j)/B) using the 
expression above is 



VB = |bx V(^-g0)0(0)+bx V0-^|^, (2) 



where 4> is the electro-static potential, b is the unit vec- 
tor along the magnetic field B, i? is the magnetic field 
strength, and the FLR effects connected with the gyro- 
average of the potential have been neglected for the sake 
of clarity. The first term in the expression above is due 
to the gradient of the eikonal, whereas the second term 
is due to the derivative of the envelope. For micro- 
instabilities, n (X 1/p* ^ 1. The second term is one 
order smaller in the normalised Larmor radius compared 
with the first and is, therefore, neglected in the lowest 
order local limit. In this paper the term 'lowest order lo- 
cal limit' will be used to indicate the gyro-kinetic model 
that retains terms only to lowest relevant order in the 
normalised Larmor radius. In this limit the radial vari- 
ation of plasma and geometry parameters is neglected, 
i.e. the model is local to a specific magnetic surface. The 
new physics, finite p*, terms introduced in this paper 
do not introduce a radial inhomogeneity and, therefore, 
the model remains 'local', but not 'lowest order'. The 
discussed physics effects are therefore distinct from any 
form of profile shearing, which relies on the introduction 
of a radial inhomogeneity. In the lowest order local limit 
the solution of the gyro-kinetic equation yields a poten- 
tial perturbation that is symmetric in 9 |15l - [l7| . It can be 
directly verified that the first term in Eq. ([2]) then yields 
a contribution that is symmetric in 9, whereas the sec- 
ond term is anti-symmetric. Based on this model we are 
searching for higher order terms that breaks the parity 
symmetry of gyro-kinetic equation The higher order p^, 
contribution to the ExB velocity breaks the parity sym- 
metry of the gyro-kinetic equation and leads to a finite, 
non-diffusive, flux of toroidal momentum. 

A simple estimate of the consequences of the higher 
order contribution in the ExB velocity on the mo- 
mentum flux can be obtained by considering the gyro- 
kinetic equation for singly charged ions retaining only 
the effect of the fluctuating ExB velocity in the back- 
ground Maxwell distribution, and the acceleration along 
the magnetic field due to the electro-static potential 



tion, and Fm the Maxwellian of the ions 
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where v\\ is the parallel velocity, e the elementary charge, 
771 is the ion mass, / the perturbed ion distribution func- 
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exp 



{v\\ - RBtu^/Bf + 2iiB/m 



(4) 

In the equation above uq is the particle density, Vth = 
^2T/m the thermal velocity, uj^{r) is the radial profile 
of the angular toroidal rotation frequency, Bt the toroidal 
component of the magnetic field, R the major radius, 
p = rav\/2B the magnetic moment, T the temperature, 
and v±_ the velocity perpendicular to the magnetic field. 
Building the parallel velocity moment of the (reduced) 
gyro-kinetic equation, assuming a Maxwell closure (for 
the moments see Ref. [s^), yields an equation for the 
perturbed parallel velocity {w) 
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where kg — nq/r, and b • V(/) = {l/qRQ)d(j)/d9 have been 
used. Furthermore, the gradient of the Maxwellian has 
been evaluated at the surface for which = (i.e. there 
is no Coriolis pinch contribution), and finite inverse as- 
pect ratio effects have been neglected {RBf/B « Rq, 
b « e^). 

The radial flux of toroidal momentum averaged over 
the flux surface is 



mnov^w, 



(6) 



where is the radial component of the ExB velocity, 
and finite inverse aspect ratio effects in the flux surface 
average have been neglected. Substituting the eikonal 
form of Eq. ([ij and integrating over the toroidal angle 
yields 



l-<l^d9Rc {2mnov'Ew'') , (7) 



where Re denotes the real part and the dagger the com- 
plex conjugate. From the equation of the perturbed par- 
allel velocity, using the complex frequency uj = ujji + 
one can derive 
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Then substituting the expression for the ExB velocity, 
and retaining the lowest order relevant terms only, 
one obtains 
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above predicts a stationary rotation gradient u'^q 
-{Rllvth){du:Jdr) 



To proceed the potential amplitude is split in a real and 
imaginary part 



+101. 



(10) 



Using the periodicity constraint in the poloidal angle, one 
obtains 

/d.^^ = 2i/d.^,§i. (11) 

At this point a particular choice for the potential 
must be made. If no symmetry breaking mechanisms are 
present the potential is symmetric in 6. However, the 
higher order contribution to the ExB velocity breaks 
the symmetry and, since the ExB velocity directly con- 
tributes to the drive of the instability, we assume a po- 
tential that has a form similar to the radial component 
of the ExB velocity 



1 Cfflo 
nq do ' 



(12) 



where 00 is symmetric in 9. We note that a finite mo- 
mentum flux is found from the equations given above 
even if we choose = 0o ■ This is because finite terms 
have been kept in the equation for the perturbed veloc- 
ity w (Eq. ([5])). The second term on the right hand side 
of Eq. ([T^ that represents the broken symmetry for the 
potential and is smaller than the first term, further 
enhances the momentum flux. The asymmetry in the 
potential, proposed in the equation above, is directly ob- 
served in the numerical simulations that will be discussed 
in Section III. 

Substituting the expression for the potential given by 
Eq. (1121) . and using the symmetry of 0o, again considering 
the lowest order relevant p* cx 1/n terms, one obtains 
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Using this expression in the equation for the momentum 
flux we finally obtain 
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where the Larmor radius is p = mvth/^B. If no ex- 
ternal torque is applied on the plasma the equation 
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The sign of the gradient means that the plasma centre 
would rotate in the counter-current direction, if the ro- 
tation frequency at the edge is zero. The ratio of the 
integrals can be estimated by introducing a scale length 
Lg for the potential perturbation 
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The effect is obviously a finite p* effect, but for standard 
parameters q = 1.4, kgp ^ 0.3, e = 0.16, and Lg — I one 
obtains 



''EQ 
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(17) 



The estimate above indicates that for the DIII-D cyclone 
base case [S^l \u'^q\ — 0.38 even without any external 
momentum input. This value is to be compared with 
the range \u'^q\ = [0 — 1] obtained in H-mode HI] with 
neutral beam heating and, consequently, with an applied 
external torque. An alternative method to estimate its 
relevance is to compare this contribution with the more 
familiar ExB shear. Adopting the estimate of the ExB 
shear rate = p*{R/LtY (for details see Ref. [11]) one 
obtains for the cyclone case \u'^q\ = 47p*. The mecha- 
nism discussed here can dominate over the ExB shearing 
but, of course, the relative strength depends on the actual 
plasma parameters considered. Even though the momen- 
tum flux introduced in this section is due to a finite p* 
correction, it is far from negligible. When developing the 
p* expansion, all quantities like e and kgp are taken to 
be of order unity, but the products of such quantities can 
still result in a factor comparable to 1 /p* . 

The mechanism discussed above is only one of several 
obtained when considering derivatives towards the par- 
allel coordinate. Following the discussion above we can 
identify several finite p* effects that might lead to sub- 
stantial momentum flux contributions: 

A Higher order p* corrections to the ExB velocity. (The 
effects of this mechanism have been estimated above.) 

B Higher order p* corrections in the treatment of the 
drift due to the magnetic field inhomogeneity (v^i) 

C Higher order p* corrections to the particle trapping. 

D Higher order p* corrections in the gyro-average as well 
as polarisation. (Note that these two are interlinked.) 

E Higher order p* correction to the calculation of the 
fluxes. (This is linked with the first point) 

It will be shown below that each of these effects breaks 
the parity symmetry of the gyro-kinetic equation, and 
leads to momentum transport. 
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II. THE MODEL 

Starting point of the model used in this paper is the 
gyro-kinetic equations of motion 



dX 
"dT 

dt 



dX 



dt 



(ZeV(0) + fiWB) 



(18) 



where is the drift velocity due to the inhomogeneous 
magnetic field, and is the ExB drift, 



vb = 



B X VV' 




b X V(0) 
B 



(19) 



In this paper we do not consider the drift due to the 
Coriolis and Centrifugal forces [H, [sl] . In Eq. (|T8l) the 
angle brackets around the potential (ip) indicate the gyro- 
average 



{4>) = ^^da</'(R + p), 



(20) 



where a is the gyro- angle, p is the Larmor radius, and 

p = p(ei cos a -f 62 sin a) (21) 

is the vector pointing from the gyro centre to the parti- 
cle position. The vectors ei and 62 are orthogonal unity 
vectors perpendicular to the magnetic field (B). The 
equations above are not correct to all orders in the nor- 
malised Larmor radius [37, 38]. The model, therefore, 
is unable to retrieve all finite corrections. It is, how- 
ever, able to describe the finite effects mentioned in 
the introduction as it will be shown below. 

Field aligned Hamada coordinates {tp, (, s) are used. 
For these coordinates the contra-variant components of 
the magnetic field are flux functions and B"^ — B'^ — 0. 
The radial coordinate ?A is a flux label, and the coordi- 
nates are chosen such that ^ remains an ignorable coor- 
dinate, i.e. all scalars that satisfy toroidal symmetry are 
not a function of C (For instance B — B{ip, s)). The co- 
ordinate s acts as a parallel coordinate B • V = B^d/ds. 
The perturbed distribution function as well as the per- 
turbed fields are assumed to have a scale length of the 
order p in the ip and C direction, but a scale length Rq 
along the magnetic field. This leads to the ordering 



O 



dl 

ds 



P*0 



P' 



dc 



(22) 



It must be noted here that, although s acts as a parallel 
coordinate, Vs does not point along the magnetic field. 
In fact, in the simplified s — a geometry s — 0/2tt and 



Vs points in the poloidal direction. The higher order 
derivatives towards s that are neglected in the lowest or- 
der local limit, therefore, are similar to the derivatives of 
the envelope in the ballooning transform. Below we will 
go through the various effects mentioned in the intro- 
duction and discuss what changes to the local model are 
necessary to retain them. The momentum fluxes are cal- 
culated using the gyro- kinetic code GKW, and for details 
on the local model the reader is referred to Ref. ISal. 



A. Higher order p« corrections to the ExB velocity 

In the linear theory the convection due to the ExB 
velocity is kept only for the background Maxwell distri- 
bution (Fa/) 



(23) 



where x = "tp, x =C, x — s, and the Einstein summa- 
tion convention has been applied. Due to the ordering, in 
the lowest order local limit, the derivative of (0) towards 
s is neglected as it results in a term that is smaller by one 
order in the normalised Larmor radius. The additional 
term, not considered in the lowest order local limit, that 
we have to add to the model in order to describe the 
flnite p* effects therefore is 



— • (Vs X W' 
B 



ds dtp 



(24) 



Here, and below, we have used the symbol — to indicate 
the finite p* terms, that have to be considered addition- 
ally to those of the lowest order local limit. Because the 
cross product of the gradients of the coordinates often 
appears in the equations, it is useful to define the tensor 



£ 



af3 



2B 



(Vx" X Wx^), 



(25) 



where the factor 2 has been introduced to make the def- 
inition equivalent with Ref. 35]. Then 



d{cP) dFM 
ds d-tp 



(26) 



B. Higher order p, correction in the drift due to 
the magnetic field inhomogeneity 

A similar argument applies to the convection caused 
by the drift velocity 



VZ3 • V = 



d 

ax"' 



(27) 



Because the derivative of any perturbed quantity towards 
s is one order smaller in the normalised Larmor radius, 
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such derivatives are neglected in the lowest order local 
limit. Therefore the additional finite term that has to 
be considered is 



VD • V 



V 



d 




OS 



C. Higher order p, corrections to the particle 
trapping 

In developing the lowest order local limit only the low- 
est order relevant p* terms are retained. The particle 
trapping is then evaluated considering the parallel con- 
vection only. For a finite beta plasma there is, however, 
a finite p* correction to the trapping connected with the 
drift motion 



a/ 

dt 



Ze 



ds dv\\ 



(29) 



This correction is not only small in p», it is also propor- 
tional to /?' and, therefore, negligible for low /? experi- 
ments. 



Because of the ordering the second term in the square 
brackets of the equation above is one order smaller in p, . 

All simulations that retain the finite p» terms in the 
gyro-average, are performed using finite difference in the 
radial and s— direction while a spectral representation is 
used for the (^-direction. The integral of the gyro-average 
is performed using 32 points on the gyro-ring. The val- 
ues of the function in between the radial grid points is 
obtained by linear interpolation. The derivative towards 
s in the equation above is calculated using central differ- 
encing on the s-grid. 

Modifications in the gyro-average directly affect the 
gyro-kinetic Poison equation: 



T 



d\{m)-c^)FM 



Ze / d^v (/) - no 



e(j) 



(33) 



where Z is the ion charge number, Te the electron tem- 
perature, and the last term on the right hand side rep- 
resents the adiabatic electron response. In this equa- 
tion the gyro-average of both the perturbed ion distribu- 
tion function / as well as the gyro-average of the gyro- 
averaged potential appears. Consistency demands that 
when the gyro-average of the potential in the evolution 
equation is modified, the same modification is applied to 
the Poisson equation. 



Higher order p* corrections in the gyro-average 
as well as polarisation 



D. 



The gyro- average is defined through Eq. (|20|) . with the 
vector p given by Eq. (PT|) . The unit vectors in the latter 
equation can be taken to be 



ei 



62 



b X W 



(30) 



V\W\ V\W\ 

The variation of the potential in the s— direction is small 



dip 



(31) 



and the gyro-average in the lowest order local limit is 
performed considering the dependence of ^ on and C 
only. In principle the full gyro-average can be performed 
by integrating over the ring retaining also the dependence 
of (j) on s. However, since the variation with s is small 
it is easier to use a Taylor expansion of (p in s, retaining 
only terms up to the first order 



(0)(V',C,s) 



1 

2^ 



da 



0(V' + p-v^,C + p- vc,s) 



ds 



p-Vs 



(32) 



E. Higher order p, terms in calculating the fluxes. 

The quasi-linear toroidal momentum flux is evaluated 

mv\\RBt 



as 



B 



■f 



(34) 



where the brackets {} denote the fiux surface average. 
Consistent with the finite p, correction of v^; • V Fm, a 
correction to toroidal momentum fiux appears 



r^;±{2/d3v£^^M!!^;|, 



(35) 



with similar expression for the fluxes of particles and en- 
ergy. 



F. Model set of equations 

For the sake of completeness, and to document exactly 
which equations are being solved, in this section the full 
set of equations is given. The notation, and normalisa- 
tion are given in Ref. [3^1, and the reader is referred to 
this paper for further details. The evolution equation for 
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the perturbed distribution / consists of several contribu- 
tions 



^ = I + II + III + IV + V + VI 
ot 

The various terms in this equation are 

1 = -Vnh-Vf 



df 

ifc„ 



II = - VZ5 • V/ ^ -^TBEDVPf+ 

Z OS 



III 



Bi 

1 df 
(t;||b + vz5) -^VB-^ ^ 

mti|| ov\\ 

df T ,,dB /3' df 



IV 



OS 



(36) 
(37) 

(38) 

(39) 

(40) 



c 
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1 2v» RBt ; 

Ln Lt Vr B 

1 1 2w|| RBt , 

Ln Lt Vr B 



F, 



M 



Fm, 



(41) 



V = 



VI = 



Ze Z B' d{(t>) 



(42) 



- -^^D ■ y{4>)FM ^ -iEoVPkp (0) Fm- 



P^EdV'^Fm, (43) 



ds 



B2 



where the tensor T> is related to £ through 



and Et is 



Ej 



-2£ 



1 dB 
Bd^ 



'II 



- 



(44) 



(45) 



Because we do not consider the non-linearity, or the neo- 
classical transport, the numbering of the terms is dif- 
ferent from Ref. [s^- Also compared with Ref. [3^1 we 
consider only the electro-static case. The Latin index p 
indicates a summation over p = 1, 2 only, i.e. excluding 
the parallel direction. The higher order corrections 
introduced in this paper are marked with under-braces. 
Their labels correspond to the enumeration of the sym- 
metry breaking mechanism in the introduction, if the 



drift velocity in the perturbed (Bi) and background {B2) 
distribution is considered as the mechanism B. 

Since is a small parameter, one can expect that 
some linearisation of the solution around = is pos- 
sible. This implies that the momentum flux generated 
by the newly introduced terms in the Eqs. pTMSp is lin- 
ear in . It also suggests that the total momentum flux 
generated by the finite p* corrections is the sum of the 
momentum fluxes generated by each of the p* terms indi- 
vidually, since their interaction would scale as p^. Both 
these hypotheses will be tested in the next section. 



G. Symmetry breaking 

The parity symmetry discussed in [l^, [l^, [l^ 
linear case considered here the transformation 



For the 



(46) 



leaves the gyro-kinetic equation invariant. Note that 
the transformation above implies for the wave vectors: 
kip — >■ — fc^, — -f fc|j. Assuming an up-down symmetric 
equilibrium, the tensors £ and V can be shown to have 
the following properties 



£^^{s) = +£'^^{-s) 



V'{s) = +V\-s) 



(47) 
(48) 
(49) 



It can then be directly verified that none of the terms re- 
tained in the lowest order local limit changes sign under 
the transformation given above. With only these terms 
kept in the model equations, the solution for the poten- 
tial is symmetric in the low field side position, while the 
parallel velocity perturbation is anti-symmetric. The re- 
sulting momentum flux is then zero. It can however, also 
be verified that all finite p* terms do change sign. All 
these terms will break the symmetry and can generate a 
finite momentum flux. 



III. RESULTS 

All simulation results in this paper are obtained with 
the gyro-kinetic code GKW 'sE], using a magnetic equi- 
librium of concentric circular surfaces where the full de- 
pendence on inverse aspect ratio e is kept [1^ , rather than 
the s—a geometry where only the lowest order in e is kept. 
For this equilibrium, j3' = 0, and the correction to the 
trapping (mechanism C in the previous section) will not 
be considered. All of the simulations have been obtained 
using finite difference in the radial and s— direction, and 
a spectral representation for the ^-direction. Periodic 
boundary conditions are used in the radial direction, con- 
sistent with the homogeneous nature of the model. The 
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radial grid, in contrast to the ballooning representation 
discussed in the introduction, corresponds to a range of 
radial wave vectors rather than one mode with a radial 
wave vector that is zero at the low field side of the torus. 

Without any of the terms in Eqs. (|37II43I) the po- 
tential {(f)), density and temperature perturbations are 
symmetric in the low field side position s — 0, whereas 
the parallel velocity perturbation (w) is anti-symmetric. 
When the finite terms are introduced the symmetry 
is broken for all quantities. The effect is clearly visi- 
ble on, for instance, (f>{s) but can nevertheless be better 
displayed by constructing the symmetric {G^) and anti- 
symmetric {G^) component of every perturbed quantity 



[Gis] + G(-s)] G' 



= i[G(s)-G(-.)] 



(50) 

Fig. [T] shows the real and imaginary part of the potential 
{(p) and parallel velocity perturbation (w) decomposed 
in their symmetric and antisymmetric parts. Shown is 
the eigenfunction averaged over the radial domain as a 
function of s ~ ^/Stt, calculated considering all finite p^, 
terms discussed in the previous section. The functions in 
the lower figure, (j)'^ and , are zero in the absence of 
finite terms. It can be seen that the finite terms 
lead to a modest asymmetry in the potential and parallel 
velocity perturbations. 



FIG. 1: (Color online) Re[(^] [solid (blue) line with symbol 
'-|-'], Im[(/)] [dashed (blue) line with symbol 'x'], Re[?j;] [solid 
(red) line with symbol 'o'], and Im[w] [dashed (red) line with 
symbol '*']. The top figure gives 0^ and w^, while the bottom 
figure gives 0"*, . All functions shown in the lower figure 
are zero in the absence of finite p, terms. 



chosen parameters are those of the Waltz standard case 
[40| : ion temperature gradient length R/Lt — 9, density 
gradient length R/Lpf = 3, magnetic shear s = 1, safety 
factor q — 2^ inverse aspect ration e = 0.16, ion to elec- 
tron temperature ratio Ti/T^ — 1. The adiabatic electron 
approximation is employed, and we consider the linear 
stability and quasi-linear fluxes only. For the chosen pa- 
rameters, the most unstable mode is the ion temperature 
gradient (ITG) mode. Unless otherwise specified, the 
normalised Larmor radius p» is 2.5 • 10"'^, and poloidal 
wave vector kep is 0.43. Note that the definition of p» is 
different by a factor \/2 compared with most literature 
due to a factor 2 the in thermal velocity vth = \/2T jm^ 
i.e. p = \plps with ps — ^JmiTe/eB. Furthermore, in p* 
the Larmor radius is normalised with the major radius 
(i?o) rather than the minor radius of the last closed flux 
surface (a). For a/R = 1/3, a/ps = 188 which is close 
to the value of the cyclone base case (s^l- Simulations 
use 30 points in the s— direction, 16 p and 16 v\\ , grid 
points. Finally, 41 radial points with a box size of 20 p 
are used in the radial direction. The representation on a 
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FIG. 2: (Color online) Re[<^"*] calculated directly from the 
eigenfunction [solid (blue) line with symbol '-I-'], Im[(/!>^] calcu- 
lated directly from the eigenfunction [dashed (blue) line with 
symbol 'x'], Re[<^"*] calculated through numerical differentia- 
tion of (f)^ [solid (green) line with symbol 'A'], and lm[0^] 
calculated through numerical differentiation of (f>^ [dashed 
(green) line with symbol 'V'] 

The eigenmode, together with the decomposition of 
Eqs. (j50p . can be used to verify the assumption repre- 
sented by Eq. (|12l) of the analytic model. This assump- 
tion gives the following relation between the symmetric 
and anti-symmetric part of the potential 

0^ = -^. (51) 
nq aO 
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Fig. [2] shows the anti-symmetric part of the potential 
calculated directly from the eigenfunction, as well as 
through numerical differentiation of the symmetric part. 
Here, only the finite correction in the ExB velocity 
is retained in the numerical simulation (i.e. only contri- 
bution A of Eqs. pTMSp is kept), in agreement with the 
assumptions made in the analytic model. Comparing the 
curves it can be seen that the agreement is reasonable, 
though not perfect. This verifies that the analytic model 
gives a reasonable estimate of the effect. Numerical sim- 
ulation is nevertheless necessary if accuracy is required. 

In this paper we consider two contributions to the mo- 
mentum fiux only: the diagonal (diffusive) contribution 
and the fiux due to the finite terms. As discussed in 
the previous section, due to the smallness of p», one can 
expect the various mechanisms to generate a momentum 
fiux that is linear in . The equation for the momentum 
fiux may then be written in the form 

= x^W + Cp, (52) 

where u' = —R^yuj^/vth- The first term in the equation 
above is the momentum diffusivity which has been stud- 
ied in several papers [H, [H, |4l| . In the absence of an 
external torque, the momentum fiux due to the finite 
parallel derivatives will modify the rotation profile, until 
this flux is balanced by the diffusive momentum flux, and 
= 0. The equilibrium toroidal rotation gradient then 

is 

(J 

u'eq = P* (53) 

The equation above provides a useful alternative to ex- 
press the magnitude of the flnite momentum fiux, as it 
more clearly indicates its impact on the rotation profile. 

Not only the fiux due to the finite terms can be 
expected to be linear in p*, the smallness of p* also sug- 
gests that the various effects are additive. Both these 
properties are shown in Fig. [31 which gives u'^q as a 
function of p* for each of the mechanisms discussed in 
the previous section, as well as the sum of all these 
separate contributions and the momentum flux calcu- 
lated including all mechanisms. It can be seen that the 
flux due to all mechanisms is nearly perfect linear in 
p*. For p* = 0.005(0.009) the momentum flux calcu- 
lated including all mechanisms is 1%(5%) smaller than 
the sum of the contributions calculated for each of the 
mechanisms separately. The various symmetry breaking 
mechanisms therefore do not have a strong interaction. 
The finite p, terms only have a small infiuence on the 
growth rate of the instability and the quasi-linear ion 
heat flux. For p, = 0.005(0.009) the growth rate in- 
creases by 0.4%(1.4%) and the ion heat fiux decreases by 
0.9%(3%). 

It can be seen from the figure that the largest con- 
tribution to the momentum transport is the fluctuating 
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FIG. 3: (Color online) The momentum flux expressed in the 
equilibrium u'^q as a function of p, . The letters on the 
right relate to the symmetry breaking mechanism identified 
in Eq. P7I43|I A: red line with symbol '-I-', Bi: light blue line 
with symbol B2: dark green with symbol 'o', D: black line 
with symbol 'V', E: magenta line with symbol 'A'. Addition- 
ally the sum of all individual contributions is given by the 
(green) curve with the symbol 'o', and the momentum flux 
calculated including all flnite p* terms is given by the (blue) 
curve with symbol 'x'. The p, axis contains typical values of 
three machines: ITER, AUG and TCV . 



ExB velocity in the background gradient, i.e. the mech- 
anism that was investigated in the introduction through 
an analytic model. Furthermore, C « 100 for the re- 
sults shown in Fig. [1] in reasonable agreement (but larger 
than) the analytic estimate C — 3/ {q{kgp)'^eLg) — 51 
(using Le ~ 1). The effect of the drift in the perturbed 
distribution on the momentum flux, as well as the cor- 
rections to the gyro-average and the calculation of the 
flux are smaller, though non- negligible, compared with 
the effect of the ExB velocity. For the parameters used 
in the simulations the sum of all effects is close to the 
contribution of the ExB velocity, due to a cancellation of 
the other contributions. This is, however, a coincidence 
for this particular set of parameters. 

Fig. 21 shows the momentum transport for different 
plasma parameters (with ~ 2.5 • 10~"^). For all param- 
eters the momentum transport due to the p* terms leads 
to a negative u'^q as predicted by Eq. Assuming 
the rotation at the plasma edge is zero, a counter-current 
rotation is then generated in the core. It can be seen that 
u'eq increases with decreasing e, in agreement with the 
prediction of the analytic model. Also the increase with 
decreasing kgp predicted by this equation is found at suf- 
flciently large kgp. Below kgp = 0.3 however, the magni- 
tude of the flux decreases with decreasing hep. Because 
also Xv decreases [13 this represents a strong decrease 
in the momentum flux due to the flnite p^ terms. This 
decrease in the flux might be related to a larger exten- 
sion of the mode along the surface (i.e. a larger Lg in 
Eq. (fTB|l ). but the exact interpretation is at present un- 
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FIG. 4: Equilibrium value of u' as a function of different plasma parameters. From left to right: magnetic shear (s), safety 
factor (g), poloidal wave vector (fcp) density gradient length R/Ln and temperature gradient length R/Lt- The three curves 
that are shown in each of the graphs correspond to different values of the inverse aspect ratio: e = 0.16 solid (red) line with 
symbol 'x', e = 0.08 dash-dotted (green) line with symbol '+', and e = 0.05 dashed (blue) line with symbol 'o'. 



known. The increase of the momentum flux with safety 
factor and magnetic shear can be understood through a 
stronger localisation of the mode (smaller Lg). For the 
lowest e, no data of u'^q is shown for the highest val- 
ues of the magnetic shear, because at higher shear the 
dominant instability is the ofi-axis ITG [43] (this mode 
also sets in for kgp > 0.7). In the lowest order local limit 
two modes exist (one shifted towards positive s and one 
towards negative s) that have equal growth rates, and 
drive an equally large momentum flux in the opposite 
direction resulting in a zero momentum flux on average. 
Even a small symmetry breaking results in one of these 
modes to be more unstable than the other. The time in- 
tegration projects out only the most unstable mode that 
than drives a large momentum flux which magnitude no 
longer depends on the magnitude of the symmetry break- 
ing term. This case is interesting, but beyond the scope 
of the present paper. 

IV. DISCUSSION 

In this paper the effect of finite terms on the toroidal 
momentum transport is investigated. Essentially, higher 
order terms connected with the derivative towards 
the 'parallel' coordinate s have been investigated. In 
the ballooning representation, these are the terms con- 



nected with the derivative of the ballooning envelope, 
rather than the eikonal. The model equations studied in 
this paper are homogeneous in the radial direction. Fi- 
nite p* effects due to the radial profiles are, therefore, 
not included. 

This paper discusses only the quasi-linear theory, and 
assumes an adiabatic electron response. Clearly, further 
study is required to assess the magnitude of the effect 
under experimentally relevant conditions. Although the 
momentum is mainly carried by the ions, kinetic electrons 
as well as finite beta effects have in the past been found to 
have a significant impact on the momentum flux [43 - l46j . 
From the analytic model one might, for instance, expect 
that the result is sensitive to the extension of the mode 
along the field line. Finally, only non-linear simulations 
can accurately determine the momentum fiuxes. These 
studies are left to future work. 

H-mode plasmas have been reported to rotate mostly 
in the co-current direction [13] , but it appears that both 
CO- and counter current rotation can occur if the plasma 
edge does not rotate strongly 49]. Experimentally, a 
transition from co- to counter current rotation is then 
observed when the density (or density gradient) exceeds 
a threshold value [l^ HH, [1^. In our model the mo- 
mentum flux due to the higher order parallel derivatives 
is directed such that the plasma core will rotate more 
strongly in the counter-current direction. 
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